
%Load files. You may need to put the full directory name here. They are in    root\TEM
folderName = 'TEM\';
fileList = dir(strcat(folderName, '*.dm4'));
clear img;

[tags, img1] = dmread(strcat(folderName, fileList(1).name));
[tags, img2] = dmread(strcat(folderName, fileList(2).name));
[tags, img3] = dmread(strcat(folderName, fileList(3).name));

%% Fig. 1D, 1E
figure(1);
clf;
x = linspace(0, 13.1362, 2048);
imagesc(x, x, img3(:, :));
axis square;
colormap gray;
caxis([2,5] * 1e4);
hold all;
a = 0.825;
plot([a, a+2], 1*[1,1], 'color', 1*[1,1,1], 'linewidth', 10);

figure(2);
scale = 2048 / 13.1362;
plot(sum(img3(:, :), 2), -(1 : 2048) / scale);
ylim([-13.1362, 0]);
yticks(1.58 * [-13:1:0]);
xlim([4.5, 9] * 1e7)

%% Fig. 1F: FFT from large area image
figure(2);
clf;
x = linspace(0, 74.31, 4096);
imagesc(x, x, img1);

axis square;
daspect([1,1,1]);
colormap gray;
caxis([2,5.5] * 1e4);
xticks([]);
yticks([]);
hold all;

a = 0.75;
plot([a, a+10], 44.25*[1,1], 'color', [1,1,1], 'linewidth', 3);



% FFT of BiO layers
% diagnostics: check which lines we are using for the fft.
imgX = img1;

% figure(1);
list = [1944, 2030, 2118, 2206, 2292, 2380, 2468, 2555, 2643, 2730, 2817, 2904, 2992, 3080, 3167, 3255, 3342, 3430, 3517, 3605, 3692, 3780];
cutsBottom = zeros(length(list), 4096);
for i = 1:length(list)

    cutsBottom(i, :) = sum(img1(list(i)-1:list(i)+17, :), 1);
    imgX(list(i)-1:list(i)+17, :) = 0;
end

listB = [1934, 1848, 1760, 1673, 1585, 1499, 1410, 1323, 1235, 1148, 1061, 974, 884, 797, 710, 622, 535, 448, 360, 273, 186, 98];
cutsTop = zeros(length(listB), 4096);
for i = 1:length(listB)
    cutsTop(i, :) = sum(img1(listB(i)-14:listB(i)-5, :), 1);
    
    imgX(listB(i)-14:listB(i)-5, :) = 0;
end

    figure(10);
    imagesc(imgX);

    colormap gray;
    caxis([2,5] * 1e4);

figure(3);
clf;
colors = jet(10);
for i = [1:3]
    plot(x, cutsBottom(i, :) + i*1.75e5, 'color', colors(i, :));
    hold all;
end
xlim([0, 38])

figure(2);
clf;
%dx = 74.31 / 4096;
dx = 74.31 * 0.9662319529 / 4096;
fM = 1/(2*dx);
freqList = linspace(-fM, fM, 4096);
freqList = freqList;
for j = [1:22]
    plot(freqList, abs(fftshift(fft(cutsTop(j, :)))) + (j-1) * 7e6, 'color', [77,177,228]/255); hold all;
    plot(freqList, -abs(fftshift(fft(cutsBottom(j, :)))) - (j-1) * 7e6, 'color', [240, 72, 65] / 255);
    hold all;
end
xlim([0, 5]);
ylim(0.125*[-1.4e9, 1e9])

xlabel('k (nm^{-1})');
ylabel('FT Amplitude');

%% Fig. 1G: Fit the FFT peaks, plot mean position vs. layer number.
figure(11);

clf;
plot(abs(fftshift(fft(cutsBottom(2, :)))));

%fit the fft peaks.
L = length(list);
peakTop = struct('A', zeros(L, 2), 'mu', zeros(L, 2), 'sigma', zeros(L, 2), 'm', zeros(L, 2), 'b', zeros(L, 2));
peakBottom = struct('A', zeros(L, 2), 'mu', zeros(L, 2), 'sigma', zeros(L, 2), 'm', zeros(L, 2), 'b', zeros(L, 2));

dx = 74.31 * 0.9662319529 / 4096;
fM = 1/(2*dx);
freqList = linspace(-fM, fM, 4096);

for i = 1:length(list)
    %do fft.
    FTop(i, :) = fftshift(fft(cutsTop(i, :)));
    FBottom(i, :) = fftshift(fft(cutsBottom(i, :)));
    
    [peakTop.A(i, :), peakTop.mu(i, :), peakTop.sigma(i, :)] = FitLorentz(freqList(2306:2330), abs(FTop(i, 2306:2330)));
    [peakBottom.A(i, :), peakBottom.mu(i, :), peakBottom.sigma(i, :)] = FitLorentz(freqList(2214:2265), abs(FBottom(i, 2214:2265)));

end
%
figure(1);
clf;
meanPeak = mean([peakTop.mu(:, 1) / sqrt(2), peakBottom.mu(:, 1)]);

yyaxis left;
plot((0:21) / 2, (peakTop.mu(:, 1) / sqrt(2) - meanPeak) / meanPeak * 100, 'o', 'color', [77,177,228]/255, 'MarkerFaceColor', [77,177,228]/255);
hold all;
plot((0:21) / 2, (peakBottom.mu(:, 1) - meanPeak) / meanPeak * 100, '^', 'color', [240, 72, 65] / 255, 'MarkerFaceColor', [240, 72, 65] / 255);
xlabel('Distance to Interface (UC)');
ylabel('\Delta k / <k> (%)');
xlim([0, 11]);
ylim([-0.10, 0.30]);

yyaxis right;
plot((0:21) / 2, peakTop.sigma(:, 1), 'o', 'color', [77,177,228]/255);
plot((0:21) / 2,  peakBottom.sigma(:, 1), '^', 'color', [240, 72, 65] / 255);
ylim([-0.01, 0.02]);

figure(2);
clf;
plot(peakTop.sigma(:, 1) / sqrt(2), 'o'); hold all;
plot(peakBottom.sigma(:, 1), 'x');

xlabel('Number of Half Unit Cells from Twist Interface');
ylabel('Peak Width (nm^{-1})');

[mean(peakTop.mu(:, 1)),std(peakTop.mu(:, 1))]
[mean(peakBottom.mu(:, 1)), std(peakBottom.mu(:, 1))]

figure(3);
clf;
errorbar(peakTop.mu(:, 1) / sqrt(2), peakTop.sigma(:, 1) / sqrt(2) / 5, 'o');
hold all;

errorbar(peakBottom.mu(:, 1), peakBottom.sigma(:, 1) / 5, 'x');
